function [cr,desc,asce,div,J,K] = get_auxillary_mat(d)
% function [cr,desc,asce,div,J,K] = get_auxillary_mat(d)
% this function is used to generate several global auxillary matrix 
% all of them are stored in dense matrix, where only necessary elements are
% stored, they means:
%           name          size              means
%           asce  :  (d+1)*(d+2)/2 x 2    the 2th and 3th nonzero entry
%           desc  :  (d+1)*(d+2)/2 x 3    the positions of nonzero entry
%            cr   :  (d+1) x (d+1)        for ease access of local dof
%           J,K   :  (d+1)*(d+2)/2 x 2    the J,K indices, I = d-J-K
%           div   :  (d+1) x (d+1)        the matrix subdivide algorithm

desc = desc_pattern(d);
cr = cr_pattern(d);
asce = asce_pattern(d);
div = sub_pattern(d);
[J,K] = indices(d);  % they are all same for different degrees, I is not same.


function cr = cr_pattern(d)
% the matrix store the indices of degree d
% much alike the function asce_pattern and desc_pattern
% the local dof index is stored in the up-right triangle.
rows = d+1;
cols = rows;
cr = zeros(rows,cols);
count = 1;
for line = 1:rows
    for row = line:-1:1
        cr(row,line-row+1) = count;
        count = count + 1;
    end
end

function [J,K] = indices(d)
% This function computes the index vectors I,J,K associated 
% with our linear ordering for the Bnet coefficients of a 
% polynomial of degree d over an arbitrary triangle.
m = (d+1)*(d+2)/2;
J = zeros(m,1);
K = J;
pos = 1;
for sum = 0:d
    idx = 0:sum;
    J(pos+idx) = sum:-1:0;
    K(pos+idx) = idx;
    pos = pos + sum + 1;
end

function div = sub_pattern(d)
% the bi-subdivision algorithm, also the simplify decasteljau algorithm
% with the paremeter be 1/2;
div = zeros(d + 1);
div(1,1) = 1;
for i = 2:d+1
    div(i,1) = div(i-1,1)*1/2;
    for j = 2:i
        div(i,j) = (div(i-1,j-1)+div(i-1,j))/2;
    end
end

function pattern = asce_pattern(d)
% this is 1d version, suit for 2d C^0 finite element method.
% this function may be execute only once in whole comput session.
% the purpose is to get a large pattern matrix for later use.
% the from is from degree d to degree (d+1),the pattern is until d
% the result is stored in three columns .
m_rows = (d+2)*(d+3)/2;
m_cols = (d+1)*(d+2)/2;
pattern = spalloc(m_rows,m_cols,2*m_rows);
idx_col = [1];
for d_row = 1:d+1
    % each cell for a degree, get the row index,the column index used the
    % row index last time;
    idx_row = (idx_col(end)+1:idx_col(end)+d_row+1);
    % put the values of j and k to corresponding entrys
    for j = 1:d_row
        pattern(idx_row(j),idx_col(j)) = d_row - j + 1;
        pattern(idx_row(j)+1,idx_col(j)) = j;
    end
    idx_col = idx_row;
end

function desc = desc_pattern(d)
% this function may be execute only once in whole comput session.
% much like the function asce_pattern.
m = (d+1)*(d+2)/2;
desc = zeros(m,3);
begin = 1;
for j = 0:d
    idx = (begin:(begin+j))';
    desc(idx,1) = idx;
    desc(idx,2) = idx + j + 1;
    desc(idx,3) = idx + j + 2;
    begin = begin + j + 1;
end
